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Abstract 

The continuum (thermodynamic) limit is the most important tool for analyzing pattern formation in 
large networks of dynamical systems with nonlocal coupling. In this limit, the solutions of the initial 
value problems for evolution equations on large discrete networks are approximated by those for the 
limiting integro-differential equations posed on continuous spatial domains. While this limiting proce- 
dure plays a pivotal role in elucidating the mechanisms of some very important effects such as chimera 
states, multistability, synchronization, and the coherence-incoherence transition, its mathematical basis 
remains little understood. In this paper, we use the combination of the ideas and results of the theory 
of graph limits and techniques for nonlinear evolution equations to provide a rigorous mathematical jus- 
tification for taking the continuum limit and to extend this method to cover many complex networks, 
for which it has not been applied before. Specifically, for dynamical networks on convergent sequences 
of simple and weighted graphs, we prove convergence of solutions of the IVPs for discrete models to 
those of the limiting continuous equations. In addition, for sequences of simple graphs that converge 
to {0, 1}— valued graphons, we show that the rate of convergence depends on the fractal dimension of 
the boundary of the support of the graph limit. The results of this paper are used to study the regions 
of continuity of chimera states and the attractors of the nonlocal Kuramoto equation on certain bipartite 
graphs. 



1 Introduction 

Systems of nonlocally coupled oscillators have been a subject of intense research during the last decade due 
to the intricate new dynamics found in these models and their importance in applications. Nonlocally cou- 
pled systems arise in modeling diverse physical, biological, and technological systems. Examples include 
neuronal models |[T6l l2l [33l. models of pattern- forming systems |[26ll34l . models of the spread of disease 
l52l|, those of coupled lasers flV ], and Josephson junctions |[30l[35l . and power networks [9], to name a few. 
Compared to systems with local (diffusive) coupling ||8l, dynamical mechanisms for pattern formation in 
nonlocal systems are much less understood. 

The continuum (thermodynamic) limit proved to be a very important and versatile tool for analyzing 
pattern formation in large networks of dynamical systems with nonlocal coupling lfT5l [Tl [37ll29l[T3ll28ll27]| . 
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In this limit, the solutions of the initial value problems (IVPs) for evolution equations on large discrete 
networks are approximated by those for the limiting integro-differential equations posed on continuous 
spatial domains. This limiting procedure plays a pivotal role in elucidating the mechanisms of some very 
important effects such as chimera states lITSi m. multistability 1371 [T3^ . synchronization, and the coherence- 
incoherence transition li28l. 



To set the stage for the forthcoming analysis of the thermodynamic limit in the nonlocally coupled sys- 
tems, we first review a representative application. In ll37ll . Wiley, Strogatz, and Girvan studied a nonlocally 
coupled system of phase oscillators 

^ i+k 

(j)i=uj+- ^ sin - , (1.1) 

j=i-k 

where (pi : M"*" — )• §^ := M/27rZ, i £ [n] := {1, 2, 3, . . . , n} is interpreted as the phase of oscillator i, 
uj is the intrinsic frequency, and the sum models the interactions between oscillator i and k of its nearest 
neighbors from each side (cf. |[T4l[T5]| ). It is assumed that the oscillators are located on a ring and indexed 
by integers Z/nZ. 

It is instructive to view ( |1.1[ ) as a system of differential equations on graph G„ = E{Gn)) with 

the vertex set V{Gn) = [n] and the edge set 

E{Gn) = e [nf ■ dist(i,j) < k] where, dist(i,j) = min{|i - j|, n - \i- 

Let Wg„ : [0, 1]^ ^ {0, 1} such that 

WcAx^y) = 1 if G E{Gn) and {x,y) G [(i - l)n--\in-^) x [(j - l)n-\jn-^). 

The plot of function WG„{x,y) in Fig. [T^ provides the pixel picture of the adjacency matrix of G„. In 
Fig.[T^ and in similar plots throughout this paper, we place the origin of the unit square in the top left corner 
of the plot to emphasize the correspondence between Wg„ and the adjacency matrix of G„. As n — )• oo, 
{W^G„} converges to Wg{x, y) (see Fig.[l]3). 

In ||37l . the analysis of the attracting states of ( |1.1[ ) employs the continuum limit of \\.\) . To this end, 
let k = rn for some fixed r G (0, 1], and send n — oo. Then in the uniformly rotating frame of coordinates 
( |1.1[ | formally becomes 

d f 

— (/>(a:, t) = J Wg{x, y) sin t) - (j){x, t)) dy as n ^ oo, (1.2) 

where (/){x, t) describes the evolution of the continuum of oscillators distributed over /. Throughout this 
paper, / denotes [0, 1], the spatial domain of the continuous limits of the dynamical networks. Equation 
(1.2 1 is called the thermodynamic limit of ( |1.1[ ). From (1.2i, it is easy to find a family of steady state 
solutions 

= qx + c, g G Z, c G M, 

so-called twisted states, and study effectively their linear stability (cf. [37 1). This information is then 
used to infer about the attracting states of the discrete system ( |1.1| ) for large n » 1. 
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Following the same lines Girnyk, Maistrenko, and Hasler studied the Kuramoto model with repelling 
coupling |.13il 

i+k 

ei = oj-n-^ ^ sin {Oj - Oi) . (1.3) 

j=i-k 

Interestingly, along with the g-twisted states in this model they found attracting states with more complex 
spatial organization, multi-twisted states, that were not present in the original model ( |1.1[ ). 

In the analysis of either model the use of the continuum limit was instrumental. However, these examples 
raise the following questions. 



(A) Does the continuum model ( |1.2[ ) truly approximate the dynamics of the discrete model ( |1.1| ) for large 

finite n? If so, in what sense the solutions of the integro-differential equation approximate those of 
([13? 

(B) How big is the class of network topologies for which one can use the continuum limit? Is it restricted to 

special graphs like A;— nearest neighbor coupling on a ring? Can it be applied for instance to the small 
world networks, the original motivation for the analysis in |[37]| ? 



The function Wg shown in Fig.[Tj5 is the limit of the functions Wg„ (Fig.[T^) representing the adjacency 
matrices of the A;-nearest neighbor family of graphs The latter is an example of a convergent graph 

sequence and Wq is the corresponding graph limit [18]. This observation hints on the relevance of the 
theory of graph limits for constructing the continuum limits for dynamical networks. In the remainder of 
this paper, we use the combination of the ideas and results of the theory of graph limits ll5ll4l lT9l 12011211 . 
and techniques for nonlinear evolution equations Ellini to provide a rigorous mathematical justification for 
taking the continuum limit and to extend this method to cover many complex networks, for which it has not 
been applied before. 

This paper is organized as follows. We review the necessary background on graph limits in Section 2. In 
Section 3, we discuss the heat equation on graphs and graph limits. Here, we extend a classical linear heat 
equation on graphs to allow nonlinear diffusion. This extension covers many dynamical networks arising 
in applications including coupled oscillator models like ( |1.1[ ). In the same section, we formally define the 
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continuum limit for dynamical networks of a convergent sequence of dense (weighted) graphs. In this limit, 
the discrete diffusion operator becomes and integral operator with the kernel representing the hmit of the 
infinite family of graphs. We show that the IVP for the limiting equation is well-posed and admits a unique 
solution in C(0, T; L°°{I)). Further, in Theorem 3.2 we specify conditions on the kernel and the initial 
conditions, which guarantee that the solutions of the IVPs remain continuous in space over subdomains of 
/. This result is used to characterize the attractors of the continuum model. In particular, we apply it to study 
the regions of continuity of the chimera states and attractors of the Kuramoto equation on the half graph (see 
Section 6). The rest of the paper is focused on studying the relation between the solutions of the IVPs for 
discrete networks and those for the corresponding limiting integro-differential equations. To this end. In 
Section 4, for sequences of simple graphs converging to {0, 1}— valued graphons, we show that the rate of 
convergence depends on the fractal dimension of the boundary of support of the graph limit. This shows 
explicitly how the geometry of the graphon affects the accuracy of the continuum limit. In Section 5, we 
analyze networks on convergent weighted graph sequences. The results of this paper are illustrated with the 
discussion of the dynamics of two concrete models: the Kuramoto-Battogtokh nonlocal system generating 
chimera states [15] and the Kuramoto equation on the half and complete bipartite graphs (cf. Section 6). 
The final section. Section 7, contains concluding remarks and the outlook on future work. 



2 Graph limits 

In this section, we review several definitions and results from the theory of graph limits that we will need 
later. In our brief tour through graph limits, we mainly follow [3| and |3r|. For the full exposition of 
this powerful theory with many diverse applications, we refer an interested reader to the pioneering papers 
by Lovasz and Szegedy |fT9l l20l . and Borgs, Chayes, Lovasz, Sos, and Vesztergombi BIH; and to the 
monograph ifTSl . 



Undirected graph G = {V{G), E{G)) without loops and multiple edges is called simple. V{G) stands 
for the set of nodes and E{G) C V{G) x V{G) denotes the edge set. 

The convergence of a sequence of simple graphs G„ = {V{Gn), E{Gn)) is defined in terms of the 
homomorphism densities 

where hom(F, G„) is the number of homomorphisms (i.e., adjacency preserving maps V{F) — )• V{Gn))- 



The probabilistic interpretation of (2.1 1 is the probability of a random map h : V{F) — )• V{Gn) to be a 
homomorphism. 

Definition 2.1. /f79]|?l/ The sequence of graphs {Gn} is called convergent ift{F, Gn) is convergent for every 
simple graph 

It turns out that the limiting object can be represented by a measurable symmetric function W : ^ L 
We recall that / stands for [0, 1]. Such functions are called graphons and the set of all graphons is denoted 
by Wo. 



It 



In the theory of graph limits, convergence in Definition 2. 1 is called left-convergence. Since this is the only convergence of 



graph sequences used in this paper, we refer to the left-convergent sequences as convergent. 
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Theorem 2.2. / Ii9l/ For every convergent sequence of simple graphs there is a graphon W G Wo such that 
t{F,Gn) ^t{F,W) := [ TT W{x^,Xj)dx VF. (2.2) 

J\V(F)\ 

Moreover, for every W G Wq f/jere « a convergent sequence of graphs satisfying (|2.2|). 



In the analytical description of graph limits, the notion of cut-norm of a graphon plays a central role. 
While this norm is not used explicitly in the analysis of the dynamical networks below, it is crucial for 
understanding of the metric properties of graphons. For any integrable function W, the cut-norm is defined 
by 



|W^||n = sup 



W{x, y)dxdy 
SxT 



where Cj stands for the set of all Lebesgue measurable subsets of /. The cut-distance between two graphons 
W and U is defined by 

6n{U,W) = inf \\U-W^\\n. 

respectively. Here, $ denotes the set of measure-preserving bijections : / — )• I and W'^{x,y) := 
W {(j){x) , cl){y)) (see |l4l[T8l for a combinatorial definition of the cut-distance between graphs). A graph 
sequence is convergent if and only if it is Cauchy in the cut-distance H . 

Graph limits are the equivalence classes of graphons 

[W] = {[/ G Wo : 5n{U, W) = 0} . 

With a customary abuse of notation, we refer to both W and [W] as graphons. The pseudo-metric (5n(-, •) 
induces the metric on x = {[W] : W G Wo}. The metric space (x, fa) is compact [ 20il . 

In spite of a rather technical definition of a graphon, it admits a simple visual interpretation. Let G be a 
simple graph on n nodes, V{G) = [n]. Define 



WGix,y) 



1, i^(^,j)eE{G)^ndix,y)e[^,i) x [i^,^) 
0, otherwise. 



(2.3) 



Then [Wq] is the graphon corresponding to G. Note that [Wg] is invariant under relabeling the nodes of G 
while Wg is not. Wg provides the pixel picture of the adjacency matrix of G. This is how we represented 
the adjacency matrix of the /c— nearest-neighbor network in Fig.[T^ (see |[T8l for more examples). 

Example 2.3. 4791 The Erdd's-Renyi graphs. Consider a sequence of random graphs G{n, p) = 
{V{G{n,p)),E{G{n,p))), V{G{n,p)) = [n] such that F {{i, j) G E{G{n,p))} = p for any G 
[n]^. This is a convergent sequence. The limit is the graphon [W] corresponding to the constant function 
W{x, y) = pfor all (x, y) G I^. The pixel picture ofWG(n,p) shown in Fig. [Tj:. The constant graphon 
W captures the limiting density of connections in G{n,p) for large n. Using the law of large numbers, it is 
easy to show that \\WG{n,p) ~p||n — s- n — t- oo. Note that {M^G(n,p)} ^^t convergent in L^-norm. This 
highlights the importance of the cut-norm in the analysis of graph limits. 
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Example 2.4. 4791/ The W-random graphs. Let W G Wo and consider Gn = {[n], E{Gn)) such that 

W{{i,j)eE{Gn)] = w(-,^- 

\n n 

This is a convergent sequence. The limit is \W]. This simple construction provides a rich family of random 
graphs, whose limits are known explicitly. This example shows that for any measurable function W E Wo, 
there is a sequence of simple graphs Gn converging to \W]. 

Example 2.5. / f79l/ The half-graphs. Let Hn,n = {V{Hn,n), E{Hn,n)) be a bipartite graph on 2n nodes 
such that 

V{Hn,n) = {1, 2, . . . , n, 1', 2', . . . , n'}, E{Hn,n) = {{i, j ) G V{Hn,n) X V{Hn,n) ■ i < j}. 

The sequence {-ffn,n} converges to the graphon [H\ where H : ^ I is the characteristic function of the 
set{{x,y) : |a; - y| > 1/2}. 



3 The formulation of the problem 

3.1 The heat equation on discrete and continuous domains 

Let Gn = {V{Gn), E{Gn), W(Gn)) be a sequence of weighted graphs, where V{Gn) = [n] and E{Gn) 
are the set of nodes and edges respectively. Symmetric matrix W{Gn) '■ [n]^ — [—1, 1] to every edge 
G E{Gn) assigns a weight 

I 0, otherwise. 
If Gn is a simple graph, W{Gn) is a {0, l}-valued matrix. 
By the nonlinear heat equation on G„ we mean 

, n 

|.W(t) = A(")5:u;(;)z?(nf)-4")), (3.1) 



IS 



where (t) = (nj") (t) , uf^ (t) , . . . , u^n^ (t)) , A^^ are scaling coefficients. Function D : M - 
Lipschitz continuous 

\D{u) - D{v)\< L\u-v\yu,v eM^. (3.2) 
Throughout this paper, we will use A^*^^ = However, other scalings may also be used. 



If D{u) = u the coupling operator on the right-hand side of (3.1 1 is the graph Laplacian, and Equation 



(3.1 1 becomes the linear heat equation on Gn- The linear heat equation has many applications in combina- 
torial problems such as random walks on graphs Q, and dynamical problems, e.g., analysis of consensus 
protocols [231. In this paper, we focus on the nonlinear heat equation, which provides the framework for a 
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large class of dynamical networks. In particular, the discrete Kuramoto equation ( 1.3 1, which we discussed 
in the Introduction, is of this type. 



In the remainder of this paper, for a class of discrete problems (3.11, we will derive and justify the 

u{x, t)= [ Wix, y)D (uiy, t) - u{x, t)) dy. (3.3) 



continuum counterpart of (3.1 1 

d 



dt 



The kernel W will be specified separately for each class of problems that we consider below. 



3.2 The well-posedness of the IVP 



Before we set out to study the relation between solutions of the discrete and continuous heat equations (3.1 1 



and (3.3 1, we first address the well-posedness of the IVP for (|3.3[). 



It is convenient to interpret the solution of the IVP for (3.3 1, u{x, t), as a vector-valued map u : [0, T 



L°°(/). Throughout this paper, we will use the bold font to denote the vector- valued function u{t) corre- 
sponding to a function of two variables u{x, t). 

Theorem 3.1. Suppose D{-) is Lipschitz continuous, W G L°°(I^), and g € L°°{I). Then for any T > 0, 
there exits a unique solution of the IVP for {3.3} u G (M.; (I)) subject to the initial condition u{0) = g. 



Proof. We rewrite the IVP for ( 3.3 1 as the integral equation 

u(t) = S + W{x, y)D {u{y, t) - u{x, t)) dydt =: [Kvl]{x, t). (3.4) 

Let Mg be a metric subspace of C(0, r; L°° (/)) (where r > will be specified later) consisting of functions 



u satisfying u(0) = g. Then (3.4i is the fixed point equation for the operator K : Mg — Mg. Without loss 
of generality, let ||M^||j^oo(/2) = 1. We show below that K is a contraction for a small r > 0. 



Indeed, let 

where L is the Lipschitz constant of D{-). For any u, v G Mg we have 

\\Ku - Kv\\j^^^ = max \\Ku{t) - Kv{t)\\^o.^j) 



(3.5) 



< max ess sup 



/x[0,i] 



\Wix, y)\ \D {uiy, t) - uix, t)) - D {v{y, t) - v^x, t))\ dydt 



< tL max <^ / \u{y,t)-v{y,t)\dy + \\u{t)-v{t)\\j^o.(i) 

< 2tL max ||u(t) — v(t)|| roo^n • 
te[o,T] ^"^ 
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Thus, by (|3.5|) we have 



\Ku — Kv\\j\[ < -llu 



(3.6) 



By the Banach contraction mapping principle, there exists a unique solution of the IVP for < \3.3\ u G 
Mg C C(0, r; L°°(/)). Using u(t) as the initial condition, the local solution can be extended to [0,2r], 
and by repeating this argument to [0, T] for any T > 0. In a similar fashion, we can prove the existence and 



uniqueness of the solution of the IVP for ( 3.3 1 on [— T, 0] for any T > 0. Furthermore, since the integrand in 



(3.4 1 is continuous as a map L°°{I) — )• L°°{I), u is continuously differentiable. Thus, we have a classical 



solution of the IVP for ( |3.3| ) on the whole real axis. 

□ 



3.3 Spatial regularity 



The classical heat equation is a parabolic partial differential equation. As other equations of this class, it has 
a strong smoothening property. Regardless of the regularity of the initial data, the solution of the IVP for the 
classical heat equation is a smooth function of the space variables for all positive times. No such mechanism 
is present in the heat equation on graph limit. Below we show that the spatial regularity of solutions of the 
IVP is determined by the regularity of graphon W and initial condition u(0). 



Theorem 3.2. Let D 



be a Lipschitz continuous function and J = {a, f3) C /. Suppose W £ 



L°° {I ) has a weak derivative -§z.W{x, y) for all x £ J and for almost all y € /; and 



ess supj^gj 



d_ 
dx 



W{;y) 



< C^. 



(3.7) 



L2{J) 



Thenforany < T < oc, all t £ [0;^]) and a < a' < f3' < /3, the solution of the IVP for (3.3) satisfied 

uit)eH\j'), J' = ia',l3'), 

provided u{0) G L°°(/) n H^{J). 



Proof. Let 



ho = - minja' — a, (3 — (3'}. 



Then for < /i < /iq, the difference quotient 



u{x + h,t) — u{x, t) 
h 



is a well-defined function on Qt = J' x [0,T]. Further, for {x,t) G Qt, ^{x,t) satisfies the following 
equation 

d 



Q^Cix,t) = jW{x,y)h-^{D{u{y,t)-u{x + h,t))-D{u{y,t)-u{x,t))}dy 
+ j^D^Wix,y)Diu{y,t)-u{x + h,t))dy, 



(3.8) 



^H^{J) stands for the Sobolev space of all Lebesgue measurable functions / on an open interval J C such that / and its 
distributional derivative fx are in i'(J) 161. 
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where 



W{x + h,y) - W{x,y) 



h 



By multiplying both sides of ( 3.8 1 by ^(x, t) and integrating both sides of the resultant equation over J' with 
respect to x, we have 

If ^£{x,tfdx = I W{x,y)h-^{D{u{y,t)-u{x + h,t))-D{u{y,t)-u{x,t))}i{x,t)dxdy 

^ J J' Ot Jj'xl 

+ / D^W{x,y)D{u{y,t) -u{x + h,t))£,{x,t)dxdy 
=■■ T1+T2. 



Without loss of generality, we assume 



\W\\L^i^l2) < 1. 



Using u S C(0, T; L°°{I)), Lipschitz continuity of D{-), and the triangle inequality, we have 

^max esssup(^^y)gj2 \D{u{y,t) -u{x,t))\ < 2L||u||c(o,T;L«=(7-)) =: C2. 



(3.9) 
(3.10) 

(3.11) 



Furthermore, using Fubini's theorem, (3.7 1, and the standard results for the difference quotients (see, e.g.. 
Theorem 5.8.3 ifTOl ). we have 



II^xW^IIl2(J'x7) < esssupy^j \\D^W\\l2^j') < Csesssup^gj \\W {■ , y)\\ ^2 (^j^ < C4, 
and, likewise, 

l|e(0)||L2(j,)<C5||u(0)|Ui(j), 

for all < /i < /lo and constants C4 and C5 independent of /i e (0, /iq). 



(3.12) 
(3.13) 



Using ( 3. 10 1 and ( 3.2 1, we bound the first term on the right hand side (|3.9 



iTil < / Lax,trdxdy = L\\m\\hu^y 
Jj'xl 

For the second term, we use ( |3.11| ), ( |3.12 i, and the Cauchy-Schwarz inequality 

iTsI < C2 [ D^^Wat) dxdy<C2\\D^,Wh2^j,^j)\\m\L2iJ') 
< C2C,\m\\LHj'y 



(3.14) 



(3.15) 



By combining (3.9 1, (3.14i, and (3.15 1, we have 



d 



C2C4, 



where inequality 2||^(i)||2,2(j/-) < ||'^(i)||^2/jM + 1 was used. Using the Gronwall's inequality, we obtain 



mm 



L2(J') 



< 



< CiWu 



The last inequality yields a uniform in /i G (0, Hq] bound on the difference quotient ||l2(j/). Using 
the properties of the difference quotients (cf. Theorem 5.8.3 ifTOl ). we conclude that u(t) G H^{J') for all 
t G (0,T]. 

□ 



4 Networks on simple graphs 



As we have seen in the example of the coupled system ( |1.1| ) and its continuum limit ( 1.2 1, graph limits can 
be used to guess a candidate for the continuum limit of dynamical networks on convergent sequences of 
dense graphs. Recall that the kernel \Wg] in ( |1-2| ), the continuum limit of \\.\\ , was, in fact, the graph limit 
of the sequence of graphs that was used in ( |1.1[ ). This leads to a question whether this is true in general. In 
this section, we initiate the study of this hypothesis for the nonlinear heat equation for certain families of 
simple graphs. Networks on weighted graphs are considered in the next section. 

In this and in the following sections, we prove that the solution of the IVP for appropriately chosen 



continuous problem (3.3 1 approximates the solutions of the discrete problems (3.1 1 when n is sufficiently 
large. We prove this result for two classes of convergent graph sequences. In this section, we consider the 
case of a sequence of simple graphs converging to a {0, l}-valued graphon, and we study a more general case 
of convergent sequences of weighted graphs in the next section. We single out networks on {0, l}-valued 
graphons for two reasons. First, many coupled oscillator models fit into this framework (see, e.g., ||37l[T3l 
and |6.2| ). Second, for this class of networks we can explicitly estimate the accuracy of approximation of 
the solutions of the discrete models by those of their continuum limits in terms of the network size and the 



geometry of the graphon of the network (cf. Theorem 4.1 1. This result is important, because it reveals the 
structural properties of the graphs shaping the accuracy of the thermodynamic limit. 

After these preliminary remarks, we turn to describe the class of networks to be studied in the remainder 
of this section. Let 1^:/^— )-{0,l}bea symmetric measurable function. We denote the support of W by 



and its boundary by dW~^ . 



For convenience, we rewrite the IVP for ( |3.3| ) 

d 



dt 



u{t,x) = J^W{x,y)D{u{y,t) -u{x,t))dy, 
9{x). 



(4.1) 
(4.2) 



u{x, 0) 

Throughout this section, to simplify presentation we assume that g{x) is an interval step function. 

Next, we define a sequence of discrete problems. To this end, we fix n G N, divide / into n subintervals 



r(n) 



0, 



1 



n 



(n) 



1 2 

n n 



An) 
■ ■ ■ 1 -^n 



n 



1 



n 



-,1 



(4.3) 



For weighted graphs, one can also define convergence by extending the notion of the homomorphism density for this case (see 
1191 for details). We do not discuss this generalization here, because for the problems that we study in this paper a simpler (and 
stronger) form of convergence, convergence in L^— norm, is sufficient (see Section 5). 
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and define a sequence of simple graphs G„ = {V{Gn), E{Gn)) such that V{Gn) = [n] and 

E{G^) = {(i, j) G [nf : (if^ x nW+ ^ $}. 



The IVP for the nonlinear heat equation on a discrete counterpart of (4.1 ), is given by 



dt 



j-.ijeEiGn) 
c/> \ I e [n\. 



(4.4) 



(4.5) 



To compare the solutions of the discrete and continuous models, it is convenient to represent the discrete 
function u^"^ = u'^\ . . . , Un"'*)"'' as a step function on / as follows 



Un{x,t) = u['^\ if X G I-"'. 



(n) 



Then Un(x, t) satisfies the following IVP 

d 



dt 



Un{t,x) 

u{x, 0) 



Wn{x, y)D {uniy, t) - Unix, t)) dy, 



II 

9n{x) 



(4.6) 

(4.7) 
(4.8) 



where 



Wn{x,y) 



1, (4")x/j"))nw^+/ 

0, otherwise 



n 



and 



9n{x) = g'f^ if X G f G [n]. 

There are many ways for approximating g{x) by gn{x)- For concreteness, we assign g^-"^ the average value 
of g{x) on U: 

g{x)dx. (4.9) 



(n) -1 

■ = n 



(n) 



Theorem 4.1. Le? u awii Un denote the vector-valued functions corresponding to the solutions of ^4A\, 



(4.2), and {4.7)-{4.9) respectively. Then for any e > and all sufficiently large n 



-(l-b-e) 

|u-Un||c(0,T;L2(7)) < C'l?^ ^ , 



(4.10) 



where constant Ci is independent ofn, and 2b = dimBdW~^ is the upper box-counting dimension ofdW^ 
(cf § 3.1,fTT]). 



Proof. Denote £,n{x, t) = Un{x, t) — u{x, t). By subtracting (4.1 1 from (4.7 1, we have 



dt 



Wnix, y) {D (uniy, t) - Unix, t)) - D t) - uix, t))} dy 



+ 1^ [ Wnix, y) - Wix, y) ) D t) - u(x, t)) dy. 



(4.11) 
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Next, we multiply both sides of (4. 1 1 1 by t) and integrate over / 

^^^^n(2;,t)^(ix = ^Wn{x,y){D{un{y,t) -Un{x,t)) - D{u{y,t) -u{x,t))}^n{x,t)dxdy 



+ / [Wnix,y) -W{x,y)) D{u{y,t) -u{x,t))S,n{x,t)dxdy. 

Ilxl 



(4.12) 



Using the Lipschitz continuity of D{-), \\W\\iaa^2-^ = 1, the triangle inequality, and the Cauchy-Schwarz 
inequality, we estimate the first term on the right-hand side of (|4.12|) 



Wnix, y) {D {uniy, t) - Unix, t)) - D {u{y, t) - u{x, t))} ^„(x, t)dxdy 



Ixl 



<L \iUy,t)-Ux,t))Ux,t)\dxdy<2L\\UmmPy (4-13) 
J Ixl 



We estimate the second term on the right-hand side of (4.12i, using the Cauchy-Schwarz inequality and 



the bound on D(-) (cf. (3.11 1) 



Wn{x, y) - W{x, y) D {u{y, t) - u{x, t)) t)dxdy 



< esssup(^^j/^4)gj2x[o,r] \D {u{y,t) - u{x,t))\ 

for some constant C2 > independent of n. 
Using ( |4.13| l and ( |4.14| ), from ( |4.12| ) we have 



Wn{x, y) - W{x, y) ] ^„(x, t)dxdy 



<C2\\W -Wn\\mP)Un\\L^i^I) (4.14) 



J^WinWh^I) < m^n\\l2^i)+2C2\\W-Wn\\mP)Un\\LHl) 



< {AL + C2\\W-Wn\\L^(^l2))\\U\l2(^j^+C2\\W-Wn\\L^py (4.15) 
" II (n) (n) 

To estimate — VFn||L2(/2), consider the set of discrete cells x that covers the boundary of 
the support of W 

J{n) = G [nf : (if^ x if^) n dW+ / 0} and C{n) = \J{n)\ . 

Using one of several equivalent definitions of the upper box-counting dimension of a subset of M" , we have 

5^0 — log 

where Ns{dW~^) is the number of cells of a (5 x 5)-mesh that intersect dW~^ (see Equation (3.12)(iv) in 
ifTTl ). Thus, for any e > and all sufficiently large n, we have 

C{n) < n2(^+^). 
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Since W and Wn coincide on all cells I^"'^ x /j"-* for which ^ J(n), for any e > and all 



sufficiently large n, we have 



\W -Wn\\l2(j2) = j {W -Wnfdxdy <C{n)n-^ < n^^Ci-''-^). 



In particular, 

||T^-VKn||L2(/2) <C3, 

for some positive constant C3 independent of n. 
By plugging (4.17 1 into ( 4.15| ), we have 



(4.16) 



(4.17) 



||?n|li2(/) < C'4||^„||^2(7) + C2||T^ - WnWl'^^p), Ci = AL + C2C3, 



dt 



(4.18) 



where positive constants C2 and C4 do not depend on n. By Gronwall's inequality, from ( |4.18| ) we have 



sup _ WUmlm) < ( llg - gnlli^d) + ^M^LJ^y^ ) exp{C74T}. (4.19) 



te[o,T] 



Finally, from ( 4.9 1 it is easy to see that 



|g-gn||i2(j) = 0{n ^] 



(4.20) 



The combination of ( |4l9l ), ( pTTe] ), and ( [4!20l ) implies ( [4TT0l ). 

□ 



5 Networks on weighted graphs 

In this section, we study a more general case of the heat equation on convergent sequences of weighted 
graphs. First, we define two graph sequences generated by a given graphon W and then we prove the 



convergence of the corresponding discrete problems to the continuum limit (4.1 1 



Throughout this section, we assume that : — )• [— 1, 1] is a symm etric measurable function. Let 
denote the partition of / into n intervals, Vn = {^f"^ ) ^ £ N} (see (4.3 1) and 



12 n 

! ? * • • ? 

n n n 



The quotient of W and Vn, denoted W/Vn, is the complete graph on n nodes 

W/Vn = {[n],[n] X [n],Wn), 
such that weights {Wn)ij are obtained by averaging W over the sets in Vn 



{Wnh = n^ I W{x,y)dxdy. (5.1) 

li X Ij 
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The second sequence of weighted graphs is obtained in a way that is similar to the construction of 



ly -random graph (cf. Example 2.4) 



M{Sn,W) = {[n],[n] x [n],Wn), (W, 



n lij 



w 



1 a 

n n 



(5.2) 



In the remainder of this section, we prove convergence of the nonlinear heat equations on W/Vn and 



M.{Sn, W) to the continuum equation on the graphon W (cf. (4.1 1). Furthermore, we show that the former 
problems correspond to the discretizations of ( |4. 1| ) using the method of Galerkin and the collocation method 
respectively, thus, relating the problem of justification of the thermodynamic limit for dynamical networks 
to two well-known numerical schemes for equations of mathematical physics. 



We first consider the IVP for the heat equation on W/V-, 



dt 



n 



nS")(0) 



(") ■ ^ r 1 



(n) 

where o- is defined in (4.9 1. 



By associating the step function n„(x, t) with u^^^ (t) (see (4.6 1), we rewrite (5.3 1 and (5.4 1 as 

d 



dt 



Un{x,t) 



Un{x,0) = 

where Wn and gn are the step functions 

Wn{x,y) 

gn{x) 



Wn{x, y)D {un{y, t) - Un{x, t)) dy, 
I 

gn{x), 



Wij for {x,y)elt^ xlf, 



g\ ' , for X G I\ ' . 



(5.3) 
(5.4) 



(5.5) 
(5.6) 



While we do not use this fact explicitly, it is instructive to note that (5.3 1 and (5.4) can be viewed as the 
Galerkin approximation of the IVP (4.1 ) and (4.2 1. Indeed, let Kn denote a finite-dimensional subspace of 

Hn = span{(pi,(p2, ■ ■ ■,<j)n}, 
where (pi = Xj(") is the characteristic function of I^'^^ = [{i — l)n~^ ,in~^). 



Replacing u{x, t) in (4. 1 1 with 



k=l 



and projecting the resultant equation on Hn, we arrive at (5.3 ). 
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Theorem 5.1. Let u and Un be the solutions of (4.1 \, (4.2), and (5.5), (5.6), respectively. Suppose W S 
L°°(/2) Wg e L°°(/). Then 



u — u 



n||c(o,r;L2(/)) Q as n oo. 



(5.7) 



Proof. By following the lines of the proof of Theorem 4. 1 (see (4.19l), for ^^(x, t) = Un{x,t) — u{x,t) 
we obtain 



2 C'lllW^ — l^n||L2(/2) 
S^V Un{t)\\L\I) < [ l|g-gn|lL2(j) + 



te[o,T] 



Co 



exp{C2T}, 



(5.8) 



where positive constants Ci and C2 are independent of n. By the Lebesgue differentiation theorem, 

Wn —J" and g„ — )- g, as n — )■ 00, 



almost everywhere on and I respectively. Thus, the statement of the theorem follows from (5.8 1. 

□ 

The heat equation on ]HI(X„, W) is analyzed in complete analogy to the IVP for W/Vn- The IVP in this 



case remains (5.5 1 and (|5.6|) modulo the definition of the step function 



W, 



n{x,y) = Wijiov ix,y) e I'/'' X I'. 



(5.9) 



We assume that y) is a bounded symmetric measurable function that is almost everywhere continuous 
on I^. Then using the observation in Lemma 2.5 131, 

Wn{x, y) — ;> W{x, y), as n — )• 00 

at every point of continuity of W, i.e., almost everywhere. Thus, by the dominated convergence theorem, 
we have 

\\W — WnWi^i^p) — ^ as n — !• 00. 

With this observation, the proof of Theorem |5 . 1 1 applies to the situation at hand. Thus, we have the following 
theorem. 



Theorem 5.2. Let u and Un be the solutions of (4.1 ), ^4.2\ , and (5.5), (5.9), (5.6), respectively. Suppose 
W G L°°(/^), g G L°°[I), and W is continuous almost everywhere on I^. Then 



|u - u„||c(o,T;L2(/)) ^ Q as 00. 



(5.10) 



6 Examples 



to 



In this section, we illustrate the results of this paper with three examples. First, we apply Theorem 3.2 
explain the regions of continuity in chimera states. Next, we discuss the attractors of the system of Kuramoto 
oscillators on the bipartite complete graph and on the half-graph. 
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Figure 2: a) The initial conditions (6.2i for the chimera state shown in b). b) A snapshot of the chimera 
state generated by (|6.1|). 



6.1 Regions of continuity of chimera states 

Chimera states are persistent patterns of coexisting regions of spatially coherent and chaotic behaviors (see 
Fig. [2^). They were found by Kuramoto and Battogtokh in the following continuum limit of a system of 
coupled phase oscillators 1.15.1 



d 

4>{x,t)=(jO+ I G{x - y) sin {(j){y,t) - cl){x,t) + a) dy. (6.1) 







dt 

Function </> : [0, 1] x M+ — := M/27rZ describes the evolution of the phase of oscillator at x G [0, 1 



The exponential kernel = exp{—K|a;|} provides nonlocal coupling between oscillators. Equation ( 6. 1 1 
was obtained using the phase reduction from the Ginzburg-Landau equation, which describes collective 
dynamics of nonlocally coupled limit cycle oscillators (cf. |15|). The sequences of discrete problems 
converging to (6.1 1 can be obtained using on of the schemes of Section[5] 



The numerical generation of the chimera state requires a careful setup, which we review next. To trigger 
a chimera state one has to start with the appropriate initial conditions, otherwise oscillators end up evolving 
in phase. Abrams and Strogatz reported that they were unable to generate chimera states from smooth initial 
conditions [IJ. Instead, one has to initialize the system with the initial condition that combines the regions 
of coherent and incoherent spatial profiles. The following initial condition was suggested by Kuramoto (cf. 
ID): 

(l){xi,0) = h{xi)ri, where = 6exp |— 30 (xi — (1/2))^| , Xi = in~^, i £ [n], (6.2) 

and ri are independent random variables drawn from the uniform distribution on (—1/2, 1/2) (see Fig.|2^). 
The values of the other parameters are k = 4, a = 1.457 (cf. OJ). Numerical integration of (|6.1[) and 



(6.2 1 with these parameter values yields persistent patterns with coexisting regions of spatially coherent and 



chaotic dynamics. A representative snapshot is shown in Fig.Gp. 



Theorem 3.2 explains the role of the initial conditions in generating chimera states. Note that function 
h{x) in ( 6.2 1 is rapidly decaying to outside a neighborhood of 1 /2 . Therefore, the initial conditions in the 
intervals Ji = (0, 0.2) and J2 = (0.8, 1) near the endpoints of the interval [0, 1] for all practical purposes 
can be viewed if they were produced by discretization of a function that is smooth over Ji and J2 (see 
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Figure 3: The graphons representing the limits of the families of the complete bipartite graphs K„ 
the half graphs Hn,n (b). 



(a) and 



Fig. For such initial conditions, Theorem 3.2 implies that the solution (/)(x, t) will remain continuous 
on Ji and J2, because H^{Ji^2) C C(Ji^2) by the Sobolev Embedding Theorem. This explains why the 
spatial profile remains coherent over Ji and J2 for positive times (see Fig. [2^). Theorem |3.2| also implies 
that it is impossible to generate chimera states starting from smooth initial data, because for such data the 
solution of the continuum limit remains continuous over the entire domain for all t > 0. This rules out 
regions of chaotic behavior in large networks, because their solutions remain close to that of the continuous 



system by Theorem 5.1 or Theorem 5.2 This explains why one can not produce chimera states from smooth 
initial conditions in [ 1 1. 



6.2 The Kuramoto equation on bipartite graphs 

To illustrate our results for networks on simple graphs, in this subsection we discuss the Kuramoto equation 
on the two families of bipartite graphs: the complete bipartite graphs {Hn,n} and the half-graphs {Kn,n}- 



The half-graph Hn^n was defined in Example 2.5 The complete bipartite graph A'n.n = ([2?^], E{Kn,n)), 
where 

E E{Kn^n) whenever l<i<n<j< 2n. 
Both graph sequences are convergent with the limits shown in Fig. [3] 

Consider the Kuramoto equation on Kn,n 

v^r\t) = l E sin(nf)(t)-^;")(i)), ^G[2n], (6.3) 

3: ij,i)eE{K„,„) 

where u\"'^ : M — S^. We consider two models for cr = 1 and a = —1. As shown below, the space 
homogeneous (synchronous) solution is stable for the a = I model and is unstable if cr = —1. 

Along with ( |6.3| ) we consider its continuum limit 

d 



Q^u{x, t) = a I K{x, y) sin {u{y, t) - u{x, t)) dy, (6.4) 
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Figure 4: Solutions of the IVP problem for the Kuramoto equation on the half- and bipartite complete 
graphs converge to the synchronous solution for cr = 1 (a) and to the step function for a = —1 (b). 



where graphon K G Wo is the limit of {Kn,n} (see Fig. [3^). Suppose u{x, 0) G C(I). By Theorem 3.2 

any t > 0, u{x, t) G C{I) where 

C{I) = {ue L°^{I) : for any open interval J C (0, 1/2) U (1/2, 1) u | j G C{J)}. 
Here, by u\j we denote the restriction of u to J. 



for 



We look for steady state solutions of (6.4i that belong to C{I). Setting the right hand side of (6.4i to 0, 
we obtain 



[ sin {u{y, t) -u{x,t))dy = 0, xG (0,1/2), 

J 1/2 



1/2 



sin (u(y, t) — u{x, t)) dy = 0, x G (1/2, 1). 



(6.5) 
(6.6) 



From (6.5 1 and (6.6 1, we find that the only steady state solutions from C{I) are the space homogeneous 

u'^{x) = c, for X G [0, 1], 



function 



and the step function 

^ fc, XG [0,1/2), 

' \ C2, X G [1/2,1], 

where constants c, ci, C2 G and |c2 — ci| = 7r/2. 



Next, we turn to the discrete model (6.3 1. The discrete counterparts of n''(x, t) and u^{x, t) are 

= cUn G M^" and v!' = (cil^,C2lI) G M^", 
where 1„ = (1,1,...,!)"^ G M". 



The linearization of (6.3 1 about u = u yields 



n 



(6.7) 
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Matrix L is the Laplacian of the half graph Hn^n 

where In is the n x n identity matrix and J„ = Inl^- As a graph Laplacian of an undirected connected 
graph, L is a symmetric positive semi-definite matrix with a simple eigenvalue |[T2]| . Thus, the space 
homogeneous solution is neutrally linearly stable for cr = 1 and is unstable when a = —1. 



The linearization of (6.3 1 about yields 



n 

which up to a sign coincides with < \6J} . Thus, u'^ is unstable if a = 1 and is neutrally stable for a = —1. 

Note that both families of synchronous and the step function steady state solutions, and n'', are 
invariant under the perturbations along the center direction span{l2n}- This suggests that both solutions 
persist under these perturbations, i.e., the synchronou s an d step steady states are stable for o" = 1 and 
a = —1 respectively. Note that the discrete model (6.3 1 has many other steady state solutions besides 



u'^ and u*. But the latter are the only two that approximate functions in C{I) and, therefore, only these 
solutions can be attractors of the discrete system for large n (cf. Theorem |3.2| ). This conclusion is consistent 
with the numerical simulations shown in Fig. |4] Numerical experiments show that the synchronous state is 
the attractor for the Kuramoto model with a = 1, while the step function is the attractor for the model with 
a = —I (see Fig.|4|5,c). 

The analysis of the Kuramoto equation on the family of half-graphs follows the same lines with obvious 
modifications. It is interesting to note, however, that Theorem |3]2] does not apply to the heat equation on the 
half-graphs. Nonetheless, the Kuramoto equation on the family of the half-graph exhibit the same attractors 
as in the case of the complete bipartite graphs. 



7 Conclusion 



The heat equation is one of the most fundamental equations of mathematical physics. On Euclidean domains, 
the heat operator is used to model many important phenomena involving diffusion, propagation, and pattern 
formation in diverse problems of physics and biology. On Riemannian manifolds, the heat equation has 
been a powerful tool for studying the topology of the underlying manifold [32 J . Its discrete counterpart, the 
heat equation on graphs plays an important role in the spectral graph theory Q. 

Motivated by the problems of the dynamics of large networks, in this paper we have introduced and 
studied the nonlinear heat equation on graphs. For convergent sequences of dense graphs, we derived the 
heat equation on graph limits and showed that it approximates the dynamics of sufficiently large networks. 
The heat equation on the graph limit is an integro-differential equation with the kernel of the integral op- 
erator given by the graphon, a symmetric bounded measurable function, which represents the limit of the 
underlying graph sequence. 

''The rigorous stability analysis requires a center manifold reduction. This is beyond the scope of this paper (see I25II23I for the 
analysis of closely related models). 
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We studied two classes of networks on convergent graph sequences: the family of networks on simple 
graphs that converge to a {0, l}-valued graphon and two families of networks on weighted graphs. In each 
case, we have proved convergence of solutions of the discrete problems on graphs to those of the continuous 
problems on graph limits. For the case of simple graphs, we have also estimated the rate of convergence in 
terms of the network size and the upper box counting dimension of the boundary of support of the graph 
limit. This result provides the link between the geometry of graphons and the accuracy of the thermodynamic 
limit approximation of large networks. For networks on weighted graphs considered in this paper, we also 
show that they can be interpreted as the discretizations of the continuum hmit via the method of Galerkin or 
the collocation method. 

The results of this paper are illustrated with two examples. First, we explain the regions of continuity in 
the chimera states, the counterintuitive solutions of the systems of coupled phase oscillators. Next, we study 
the attractors in the Kuramoto coupled oscillators equation model on the complete bipartite graph and on 
the half graph. Both examples highlight the qualitative differences between solutions of the classical heat 
equation and that on the graph limit. While the former are smooth functions of the space variables for all 
positive times, due to the smoothening property of the parabolic partial differential operators, the latter may 
be discontinuous, as the examples of the heat equation on the complete bipartite graph shows. 

The theory of graph limits provides a useful set of tools for studying dynamics of large networks. On one 
hand, known graph Umits for various convergent sequences Uke the half graph or VF-random graphs suggest 
continuum limits for a broad class of large networks. On the other hand, this rich theory offers many useful 
ideas and analytical results that can be applied to dynamical networks. In this paper, we analyzed three 
families of networks on convergent sequences of deterministic graphs. In the future work, we will consider 
networks on convergent sequences of random graphs such as the Erdos-Renyi and the small world graphs, 
which are ubiquitous in apphcations 1(11 1241 [36]| . 
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